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Abstract. The coronal magnetic field is an important quantity because the magnetic 
field dominates the structure of the solar corona. Unfortunately direct measurements of 
coronal magnetic fields arc usually not available. The photospheric magnetic; field is mea- 
sured routinely with vector magnetographs. These photospheric measurements are ex- 
trapolated into the solar corona. The extrapolated coronal magnetic field depends on as- 
sumptions regarding the coronal plasma, e.g. force- freeness. Force-free means that all non- 
magnetic forces like pressure gradients and gravity arc neglected. This approach is well 
OO justified in the solar corona due to the low plasma beta. One has to take care, however, 

about ambiguities, noise and non-magnetic forces in the photosphere, where the mag- 
^— ^ netic field vector is measured. Here we review different numerical methods for a nonlin- 

^ ^ ear force-free coronal magnetic field extrapolation: Grad- Rubin codes, upward intcgra- 

tion method, MHD-relaxation, optimization and the boundary element approach. We briefly 
^ discuss the main features of the different methods and concentrate mainly on recently 

developed new codes. 
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1. Introduction 

Information regarding the coronal magnetic field are im- 
portant for space weather application like the onset of 
flares and coronal mass ejections (CMEs). Unfortunately 
we usually cannot measure the coronal magnetic field di- 
rectly, although recently some progress has been made 
[see e.g.. Judge, 1998; Solanki et al, 2003; Lin et al., 
2004]. Duo to the optically thin coronal plasma direct 
measurements of the coronal magnetic field have a line- 
of-sight integrated character and to derive the accurate 
3D structure of the coronal magnetic field a vector to- 
mographic inversion is required. Corresponding feasibil- 
ity studies based on coronal Zeeman and Hanle effect 
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measurements have recently been done by Kramar et al. 
[2006] and Kramar and Inhester [2006]. These direct 
measurements are only available for a few individual cases 
and usually one has to extrapolate the coronal magnetic 
field from photospheric magnetic measurements. To do 
so, one has to make assumptions regarding the coronal 
plasma. It is helpful that the low solar corona is strongly 
dominated by the coronal magnetic field and the mag- 
netic pressure is orders of magnitudes higher than the 
plasma pressure. The quotient of plasma pressure p and 
magnetic pressure, _B^/(2/i()) is small compared to unity 
{(3 = 2^(,p/B^ <^ 1). In lowest order non-magnetic forces 
like pressure gradient and gravity can be neglected which 
leads to the force-free assumption. Force free fields are 
characterized by the equations: 



j X B = (1) 
VxB = Moj, (2) 
V-B = 0, (3) 

where B is the magnetic field, j the electric current den- 
sity and /if) the permeability of vacuum. Equation (1) 
implies that for force-free fields the current density and 
the magnetic field are parallel, i.e. 

Hoi = aB, (4) 
or by replacing j with Eq. (2) 

V X B = qB, (5) 

where a is called the force-free function. To get some 
insights in the structure of the space dependent function 
Q, we take the divergence of Eq. (4) and make use of 
Eqs. (2) and (3): 

B • Va = 0, (6) 

which tells us that the forcc-froo function a is con- 
stant on every field line, but will usually change from one 
field line to another. This generic case is called nonlinear 
force-free approach. 

Popular simplifications are a = (current free poten- 
tial fields, sec e.g., Schmidt [1964]; Semel [1967]; Schat- 
ten et al. [1969]; Sakurai [1982]) and a = constant (linear 
force-free approach, see e.g., Nakagawa and Raadu [1972]; 
Chiu and Hilton [1977]; Seehafer [1978]; AUssandrakis 
[1981]; Seehafer [1982]; Semel [1988]). These simplified 
models have been in particular popular due to their rel- 
ative mathematical simplicity and because only line-of- 
sight photospheric magnetic field measurements are re- 
quired. Linear force-free fields still contain one free global 
parameter a, which can be derived by comparing coro- 
nal images with projections of magnetic field lines (e.g., 
Carcedo et al. [2003]). It is also possible to derive an av- 
eraged value of a from transverse photospheric magnetic 
field measurements (e.g. Pevtsov et al. [1994]; Wheailand 
[1999]; Leka and Skumanich [1999]). Despite the popu- 
larity and frequent use of these simplified models in the 
past, there are several limitations in these models (see 
below) which ask for considering the more sophisticated 
nonlinear force-free approach. 

Our aim is to review recent developments of the ex- 
trapolation of nonlinear force- free fields (NLFFF). For 
earlier reviews on force-free fields we refer to [Sakurai, 
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1989; Aly, 1989; Aman et al, 1997; McClymont et al, 
1997] and chapter 5 of Aschwanden [2005]. Here we 
will concentrate mainly on new developments which took 
place after these earlier reviews. Our main emphasis is to 
study methods which extrapolate the coronal magnetic 
field from photospheric vector magnetograms. Several 
vector magnetographs are currently operating or planed 
for the nearest future, e.g., ground based: the solar 
flare tclcscope/NAOJ [Sakurai et al, 1995], the imaging 
vector magnctograph/MEES Observatory [Mickey et al., 
1996], Big Boar Solar Observatory, Infrared Polarimctcr 
VTT, SOLIS/NSO [Henney et al, 2006] and space born: 
Hinodc/SOT [Shirmzu, 2004], SDO/HMI [Borrero et al., 
2006]. Measurements from these vector magnetograms 
will provide us eventually with the magnetic field vec- 
tor on the photosphere, say B^o for the normal and B^o 
and Byo for the transverse field. Deriving these quanti- 
ties from the measurements is an involved physical pro- 
cess, which includes measurements based on the Zccman 
and Hanlc effect, the inversion of Stokes profiles (e.g., 
LaBonte et al. [1999]) and removing the 180° ambiguity 
(e.g., Metcalf [1994]; Metcalf et al. [2006]) of the hori- 
zontal magnetic field component. Special care has to be 
taken for vector magnetograph measurements which are 
not close to the solar disk, when the line-of-sight and nor- 
mal magnetic field component are far apart (e.g. Gary 
and Hagyard [1990]). For the purpose of this paper we 
do not address the observational methods and recent de- 
velopments and problems related to deriving the photo- 
spheric magnetic field vector. We rather will concentrate 
on how to use the photospheric Bxo, Byo and B^o to de- 
rive the coronal magnetic field. 

The transverse photospheric magnetic field {Bxo, Byo) 
can be used to approximate the normal electric current 
distribution by 



dByo dBxo /„x 

and from this one gets the distribution of a on the 
photosphere by 



a{x,y)=no^ (8) 

By using Eq. (8) one has to keep in mind that rather 
large uncertainties in the transverse field component and 
numerical derivations used in (7) can cumulate in signifi- 
cant errors for the current density. The problem becomes 
even more severe by using (8) to compute a in regions 
with a low normal magnetic field strength Bzo- Special 
care has to be taken at photospheric polarity inversion 
lines, i.e. lines along which B^o = [sec e.g., Cuperman 
et al, 1991]. The nonlinear forcc-frcc coronal magnetic 
field extrapolation is a boundary value problem. As we 
will see later some of the NLFFF-codes make use of (8) 
to specify the boundary conditions while other methods 
use the photospheric magnetic field vector more directly 
to extrapolate the field into the corona. 

Pure mathematical investigations of the nonlinear 
force-free equations [sec e.g. Aly, 1984; Boulmezaoud and 
Amari, 2000; Aly, 2005] and modelling approaches not 
based on vector magnetograms arc important and occa- 
sionally mentioned in this paper. A detailed review of 
these topics is well outside the scope of this paper, how- 
ever. Some of the model-approaches not based on vector 
magnetograms are occasionally used to test the nonlinear 
force-free extrapolation codes described here. 
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1.1. Why do we need nonlinear force- free fields? 

• A comparison of global potential magnetic field 
models with TRACE-imagcs by Schnjver et al, [2005] re- 
vealed that significant nonpotentially occurs regulary in 
active regions, in particular when new flux has emerged 
in or close to the regions. 

• Usually a changes in space, even inside one active 
region. This can be seen, if we try to fit for the op- 
timal linear force-free parameter a by comparing field 
lines with coronal plasma structures. An example can be 
seen in Wiegelmann and Neukirch [2002] where stereo- 
scopic reconstructed loops by Aschwanden et al. [1999] 
have been compared with a linear force-free field model. 
The optimal value of a changes even sign within the in- 
vestigated active regions, which is a contradiction to the 
a = constant linear force-free approach (see Fig. 1). 

• Photosphcric q: distributions derived from vector 
magnetic field measurements by Eq. (8) show as well 
that a usually changes within an active region [see, e.g. 
Regnier et al., 2002]. 

• Potential and linear force-free fields are too simple 
to estimate the free magnetic energy and magnetic topol- 
ogy accurately. The magnetic energy of linear force-free 
fields is unbounded in a halfspacc [Seehafer, 1978] which 
makes this approach unsuitable for energy approxima- 
tions of the coronal magnetic field. Potential fields have 
a minimum energy for an observed line-of-sight photo- 
spheric magnetic field. An estimate of the excess of en- 
ergy a configuration has above that of a potential field is 
an important quantity which might help to understand 
the onset of flares and coronal mass ejections. 

• A direct comparison of measured fields in a newly 
developed active region by Solanki et al. [2003] with ex- 
trapolations from the photosphere with a potential, lin- 
ear and nonlinear force-free model by Wiegelmann et al. 
[2005b] showed that nonlinear fields are more accurate 
than simpler models. Fig. 2 shows some selected mag- 
netic field lines for the original measured field and ex- 
trapolations from the photosphere with the help of a po- 
tential, linear and nonlinear force-free model. 

These points tell us that nonlinear force- free modelling 
is required for an accurate reconstruction of the coronal 
magnetic field. Simpler models have been used frequently 
in the past. Global potential fields provide some informa- 
tion of the coronal magnetic field structure already, e.g. 
the location of coronal holes. The generic case of force- 
free coronal magnetic field models are nonlinear force- 
free fields, however. Under generic we understand that 
a can (and usually will) change in space, but this ap- 
proach also includes the special cases a = constant and 
Q = 0. Some active regions just happen to be more po- 
tential (or linear force-free) and if this is the case they 
can be described with simpler models. Linear force-free 
models might provide a rough estimate of the true 3D 
magnetic field structure if the nonlinearity is wealc. The 
use of simpler models was often justified due to limited 
observational data, in particular if only the line-of-sight 
photospheric magnetic field has been measured. 

While the assumption of norflinear force-free fields is 
well accepted for the coronal magnetic fields in active 
regions, this is not true for the photosphere. The pho- 
tospheric plasma is a finite (3 plasma and nonmagnetic 
forces like pressure gradient and gravity cannot be ne- 
glected here. As a result electric currents have a com- 
ponent perpendicular to the magnetic field, which con- 
tradicts the force-free assumption. We will discuss later, 
how these difiiculties can be overcome. 
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2. Nonlinear force-free codes 

Different methods have been proposed to extrapolate 
nonhncar forcc-frcc fields from photospheric vector mag- 
netic field measurements. 

1. The Grad- Rubin method was proposed for fusion 
plasmas by Grad and Rubin [1958] and first applied to 
coronal magnetic fields by Sakurai [1981]. 

2. The upward integration method was proposed by 
Nakagawa [1974] and encoded by Wu et al. [1985]. 

3. The MHD relaxation method was proposed for gen- 
eral MHD-equilibria by Chodura and Schlueter [1981] and 
applied to force-free coronal magnetic fields by Mikic and 
McClymont [1994]. 

4. The optimization approach was developed by 
Wheatland et al. [2000]. 

5. The boundary element (or Greens function like) 
method was developed by Yan and Sakurai [2000]. 

2.1. Grad-Rubin method 

The Grad-Rubin method reformulates the nonlinear 
force-free equations in such a way, that one has to solve 
a well posed boundary value problem. This makes this 
approach also interesting for a mathematical investiga- 
tion of the structure of the nonlinear force- free equations. 
Bineau [1972] demonstrated that the used boundary con- 
ditions (vortical magnetic field on the photosphere and a 
distribution at one polarity) ensure, at least for small 
values of a and weak nonlincarities the existence of a 
unique nonlinear force-free solution. A detailed analysis 
of the mathematical problem of existence and uniqueness 
of nonlinear force-free fields is outside the scope of this 
review and can be found e.g. in Amari et al. [1997, 2006]. 

The method first computes a potential field, which can 
be obtained from the observed line-of-sight photospheric 
magnetic field (say Bz in Cartesian geometry) by differ- 
ent methods, e.g., a Greens function method as described 
in Aly [1989]. It is also popular to use linear force- free 
solvers, e.g. as implemented by Seehafer [1978]; Alissan- 
drakis [1981] with the linear force-free parameter a = to 
compute the initial potential field. The transverse com- 
ponent of the measured magnetic field is then used to 
compute the distribution of a on the photosphere by Eq. 
(8). While a is described this way on the entire photo- 
sphere, for both polarities, a well posed boundary value 
problem requires that the a distribution becomes only de- 
scribed for one polarity. The basic idea is to itcratively 
calculate a for a given B field from (6), then calculate 
the current via (4) and finally update B from the Biot- 
Savart problem (5). These processes are repeated until 
the full current as prescribed by the a-distribution has 
been injected into the magnetic field and the updated 
magnetic field configuration becomes stationary in the 
sense that eventually the recalculation of the magnetic 
field with Amperes law does not change the configura- 
tion anymore. To our knowledge the Grad-Rubin ap- 
proach has been first implemented by Sakurai [1981]. a 
has been prescribed on several nodal points along a num- 
ber of magnetic field lines of the initial potential field. 
The method used a finite-element-like discretization of 
current tubes associated with magnetic field lines. Each 
current tube was divided into elementary current tubes 
of cylindrical shape. The magnetic field is updated with 
Ampere's law using a superposition of the elementary 
current tubes. The method was limited by the number 
of current-carrying field lines, nodal-points and the cor- 
responding number of nonlinear equations {N^) to solve 
with the available computer resources more than a quar- 
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ter century ago. 

Computer resources have increased rapidly since the 
first NLFFF- implementation by Sakurai [1981] and about 
a decade ago the Grad-Rubin method has been im- 
plemented on a finite difiierence grid by Amari et al. 
[1997, 1999]. This approach decomposes the equations 
(1-3) into a hyperbolic part for evolving a along the mag- 
netic field lines and an elliptic one to iterate the updated 
magnetic field from Amperes law. For every iteration 
step k one has to solve iteratively for: 



a^'^lsi = ao± (10) 

(11) 

which evolves a in the volume and 

V X B^'^+i' = a^'^B^'^^ (12) 

Bi''+'^\s± = S,o (14) 

lim |b''=+''| = 0, (15) 

|r|^oo 



where ao± corresponds to the photospheric distribu- 
tion of a for either on the positive or the negative polar- 
ity. The Grad-Rubin method as described in Amari et al. 
[1997, 1999] has been applied to investigate particular 
active regions in Bleybel et al. [2002] and a comparison 
of the extrapolated field with 2D projections of plasma 
structures as seen in Ha, EUV and X-ray has been done 
in Regnier et al. [2002]; Regnier and Amari [2004]. The 
code has also been used to investigate mutual and self 
helicity in active regions by Regnier et al. [2005] and to 
fiaring active regions by Regnier and Canfield [2006]. 

A similar approach as done by Sakurai [1981] has been 
implemented by Wheatland [2004]. The implemented 
method computes the magnetic field directly on the nu- 
merical grid from Ampere's law. This is somewhat sim- 
pler and faster as Sakurai's approach which required solv- 
ing a large system of nonlinear equations for this aim. 
The implementation by Wheatland [2004] has, in partic- 
ular, been developed with the aim of parallelization. The 
parallelization approach seems to be effective due to a 
limited number of inter-process communications. This is 
possible because as the result of the linearity of Ampere's 
law the contributions of the different current carrying 
field lines are basically independent from each other. In 
the original paper Wheatland reported problems for large 
currents on the field lines. These problems have been re- 
lated to an error in current representation of the code 
and the corrected code worked significantly better, see 
also Schrijver et al. [2006] . The method has been further 
developed in Wheatland [2006] . This newest Wheatland- 
implementation scales with the number of grid points N'^ 
for a volume, rather than for the earlier Wheat- 
land [2004] implementation. The main new development 
is a faster implementation of the current-field iteration. 
To do so the magnetic field has been separated into a 
current-free and a current carrying part at each iteration 
step. Both parts arc solved using a discrete Fast Fourier 
Transformation, which imposes the required boundary 
conditions implicitly. The code has been parallelized on 
shared memory distributions with OpenMP. 

Amari et al. [2006] developed two new versions of their 
Grad-Rubin code. The first version is a finite difference 
method and the code was called 'XTRAPOL'. This code 
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prescribe the coronal magnetic field with the help of a 
vector potential A. The code has obviously it's heritage 
from the earlier implementation of Amari et al. [1999], 
but with several remarkable differences: 

• The code includes a divergence cleaning routine, 
which takes care about V-A = 0. The condition V-A = 
is fulfilled with high accuracy in the new code 10~® com- 
pared to lO"'^ in the earlier implementation. 

• The lateral and top boundaries are more flexible 
compared to the earlier implementation and allow a finite 
Bn and non zero a-values for one polarity on all bound- 
aries. This treats the whole boundary (all six faces) as a 
whole. 

• The slow current input as reported for the earlier 
implementation, which lead to a two level iteration, has 
been replaced. Now the whole current is injected at once 
and only the inner iteration loop of the earlier code re- 
mained in the new version. 

• The computation of the q characteristics has been 
improved with an adaptive Adaiiis-Basliforth integration 
scheme [see Press, 2002, for details]. 

• The fixed number of iteration loops have been re- 
placed by a quantitative convergence criterium. 

In the same paper Amari et al. [2006] introduced an- 
other Grad-Rubin approach based on finite elements, 
which they called 'FEMQ'. Different from alternative im- 
plementation this code does not use a vector potential but 
iterates the coupled divergence and curl system, which 
is solved with the help of a finite element discretization. 
The method transforms the nonlinear force- free equations 
into a global linear algebraic system. 

Inhester and Wiegelmann [2006] implemented a Grad- 
Rubin code on a finite element grid with staggered field 
components (see Yee [1966]) which uses discrete Whitney 
forms {Bossavit [1988]). Whitney forms allow to trans- 
form standard vector analysis (as the differential opera- 
tors gradient, curl and divergence) consistently into the 
discrete space used for numerical computations. Whit- 
ney forms contain four types of finite elements (form 0- 
3). They can be considered as a discrete approximation 
of differential forms. The finite element base may con- 
sist of polynominals of any order. In its simplest form, 
the 0-forms have as parameters the function values at the 
vertices of the cells and are linearly interpolated within 
each cell. 1-forms are a discrete representation of a vector 
field defined on the cell edges. 2-forms are defined as the 
field component normal to the surfaces of the cells. The 
3-forms arc finite volume elements for a scalar function 
approximation, which represents the average of a scalar 
over the entire cell. The four forms are related to each 
other by GRAD (0 to 1 form), CURL (1 to 2 form) and 
DIV (2 to 3 form). As for continous differential forms, 
double differentiation (CURL of GRAD, DIV of CURL) 
give exactly zero, independent of the numerical precision. 
A dual grid, shifted by half a grid size in each axis, was in- 
troduced in order to allow for Laplacians. Whitney forms 
on the dual grid are related to forms on the primary grid 
in a consistent way. 

The Grad-Rubin implementation uses a vector poten- 
tial representation of the magnetic field, where the vector 
potential is updated with a Poisson equation in each it- 
eration step. The Poisson equation is effectively solved 
with the help of a multigrid solver. The main computing 
time is spend to distribute a along the field lines with (6). 
This seems to be a general property of Grad-Rubin imple- 
mentations. One can estimate the scaling of (6) by oc N^, 
where the number of field lines to compute is oc and 
the length of a field line oc N. The Biot-Savart step (5) 
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solved with FFT or multigrid methods scales only with 
oc A'^ log A'^. Empirical tests show that the number of it- 
eration steps until a stationary state is reached does not 
depend on the number of grid points A'^ for Grad-Rubin 
solvers. We have explained before, that the Grad-Rubin 
implementation requires the prescription of a only for one 
polarity to have a well posed mathematical problem. The 
Inhester and Wiegelmann [2006] implementation allows 
those choice of boundary conditions as a special case. In 
general one docs not need to make the distinction be- 
tween (dV)^ and (dV)^ in the now implementation. A 
well posed mathematical problem is still ensured, how- 
ever, in the following way. Each boundary value of a is 
attached with a weight. The final version of a on each 
field line is then determined by a weighted average of the 
a-values on both endpoints of a field lines. By this way 
the influence of uncertain boundary values, e.g. on the 
side walls and imprecise photospheric measurements can 
be suppressed. 

2.2. Upward integration method 

The basic equations for the upward integration method 
(or progressive extension method) have been published 
already by Nakagawa [1974] and a corresponding code 
has been developed by Wu et al. [1985, 1990a]. The up- 
ward integration method is a straight forward approach 
to use the nonlinear force-free equations directly to ex- 
trapolate the photospheric magnetic field into the corona. 
To do so one reformulates the force- free equations (1-3) in 
order to extrapolate the measured photospheric magnetic 
field vector into the solar corona. 

As a first step the magnetic field vector on the lower 
boundary Bo (a;, y, 0) is used to compute the ^-component 
of the electric current pojzo with Eq. (7) and the pho- 
tospheric a-distribution (say ao) by Eq. (8). With the 
help of Eq. (4) we calculate the x and j/-component of 
the current density 



U'ojxo = cuo B^o, (16) 
Mojj/O = aoByo. (17) 

We now use Eq. (3) and the x and y-component of 
Eq. (2) to obtain expressions for the 2:-derivatives of all 
three magnetic field components in the form 



dBxo ■ , dBzo 

dByo dBzo 

dBzo _ dBxo dByo 
dz dx dy 



(18) 

(19) 
(20) 



The idea is to integrate this set of equations numer- 
ically upwards in z by repeating the previous steps at 
each height. As a result wc get in principle the 3D mag- 
netic field vector in the corona. While this approach is 
straight forward, easy to implement and computational 
fast (no iteration is required), a serious drawback is that 
it is unstable. Several authors [e.g., Cuperman et al, 
1990; Amari et al., 1997] pointed out that the formular 
tion of the force-free equations in this way is unstable 
because it is based on an ill-posed mathematical prob- 
lem. In particular one finds that exponential growth of 
the magnetic field with increasing height is a typical be- 
haviour. What makes this boundary value problem ill- 
posed is that the solution does not depend continuously 
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on the boundary data. Small changes or inaccuracies in 
the measured boundary data lead to a divergent extrap- 
olated field [see Low and Lou, 1990, for a more detailed 
discussion]. As pointed out by Low and Lou meaning- 
ful boundary conditions are required also on the outer 
boundaries of the computational domain. It is also pos- 
sible to prescribe open boundaries in the sense that the 
magnetic field vanishes at infinity. This causes an addi- 
tional problem for the upward integration method, be- 
cause the method transports information only from the 
photosphere upwards and docs not incorporate boundary 
information on other boundaries or at infinity. Attempts 
have been made to regularize the method [e.g., Cuperman 
et al, 1991; Demoulin and Priest, 1992], but cannot be 
considered as fully successful. 

Wu et al. [1990b] compared the Grad-Rubin method 
in the implementation of Sakurai [1981] with the up- 
ward integration method in the implementation of Wu 
et al. [1990a] ^ The comparison showed qualitatively 
similar results for extrapolations from an observed mag- 
netogram, but quantitatively differences. The NLFFF- 
computations have been very similar to potential field 
extrapolations, however, too. One reason for this be- 
haviour was, that the method of Sakurai [1981] is limited 
to small values of a and an 'by eye' comparison shows 
that the corresponding NLFFF field is very close to a 
potential field configuration. The field computed with 
the upward integration method deteriorated if the height 
of the extrapolation exceeded a typical horizontal scale 
length. 

The upward integration method has been recently re- 
examined by Song et al. [2006] who developed a new for- 
mulation of this approach. The new implementation uses 
smooth continuous functions and the equations are solved 
in asymptotic manner iteratively. The original upward 
integration equations arc reformulated into a set of ordi- 
nary differential equations and uniqueness of the solution 
seems to be guarantied at least locally. While Demoulin 
and Priest [1992] stated that 'no further improvement 
has been obtained with other types of smoothing func- 
tions' the authors of Song et al. [2006] point out that the 
transformation of the original partial differential equa- 
tions into ordinary ones eliminates the growing modes 
in the upward integration method, which have been re- 
ported before in Wu et al. [1990a] and subsequent pa- 
pers. The problem that all three components of the pho- 
tosphcric magnetic field and the pliotosphcric a distribu- 
tion has to be prescribed in a consistent way remains in 
principle, but some compatibility conditions to compute 
a slowly varying a have been provided by Song et al. 
[2006]. These compatibility conditions are slightly dif- 
ferent for real photospheric observations and tests with 
smooth boundaries extracted from semi-analytic equilib- 
ria. For the latter kind of problems the new formulation 
provided reasonable results with the standard test equi- 
librium found by Low and Lou [1990]. The method seems 
to be also reasonable fast. Of course further tests with 
more sophisticated equilibria and real data are necessary 
to evaluate this approach in more detail. 

2.3. MHD relaxation 

MHD relaxation codes [e.g., Chodura and Schlueter, 
1981] can be applied to solve nonlinear force-free fields as 
well. The idea is to start with a suitable magnetic field 
which is not in equilibrium and to relax it into a force- 
free state. This is done by using the MHD equations in 
the following form: 
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i/v = (V X B) X B (21) 
E + V X B = (22) 

-V X E (23) 



9B 

dt 

V-B = 0, (24) 



where v is a viscosity and E the electric field. As 
the MHD-releixation aims for a quasi physical temporal 
evolution of the magnetic field from a non-equilibrium to- 
wards a (nonlinear forcc-frcc) equilibrium this method is 
also called 'evolutionary method' or 'magncto-frictional 
method'. The basic idea is that the velocity field in the 
equation of motion (21) is reduced during the relaxation 
process. Ideal Ohm's law (22) ensures that the magnetic 
connectivity remains unchanged during the relaxation. 
The artificial viscosity u plays the role of a relaxation 
coeflScient which can be chosen in such way that it accel- 
erates the approach to the equilibrium state. A typical 
choice is 



u = i |B|^ (25) 

with /i = constant. Combining Eqs. (21), (22), (23) 
and (25) we get an equation for the evolution of the mag- 
netic field during the relaxation process: 



c?B 

= M Fmhd (26) 



with 



F.Ho = Vx( t(^x^)^:^]x^ ). (27) 

This equation is then solved numerically starting with 
a given initial condition for B, usually a potential field. 
Equation (26) ensures that Eq. (24) is satisfied during the 
relaxation if the initial magnetic field satisfies it. ^ The 
difficulty with this method is that it cannot be guaran- 
teed that for given boundary conditions and initial mag- 
netic field (i.e. given connectivity), a smooth force-free 
equilibrium exists to which the system can relax. If such 
a smooth equilibrium does not exist the formation of cur- 
rent sheets is to be expected which will lead to numerical 
difficulties. Therefore, care has to be taken when choos- 
ing an initial magnetic field. 

Yang et al. [1986] developed a magneto frictional 
method which represent the magnetic field with the help 
of Euler (or Clebsch) potentials. 



B = Vff X V/i, (28) 

where the potentials g and h arc scalar functions. The 
general method has been developed for three dimensional 
fields and iterative equations for g{x, y, z) and h{x, y, z) 
have been derived. The Clebsch representation automat- 
ically ensures V • B = 0. The method has been explicitly 
tested in the paper by Yang et al. [1986] with the help 
of an equilibrium with one invariant coordinate. In prin- 
ciple it should be possible to use this representation for 
the extrapolation of nonlinear force- free fields, but we are 
not aware of a corresponding implementation. Due to the 
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discussion in Yang et al. [1986] a difficulty seems to be 
that one needs to specify boundary conditions for the po- 
tentials, rather than for the magnetic fields. It seems in 
particular to be difficult to find boundaries conditions for 
potentials which correspond to the transverse component 
of the photospheric magnetic field vector. One problem 
is that boundary conditions for g and h prescribe the 
connectivity. Every field line can be labelled by its {g, h) 
values. Hence boundary values for g and h establish foot 
point relations although the field is not known yet. 

The MHD-rclaxation (or evolutionary) method has 
been implemented by Mikic and McClymont [1994]; Mc- 
Clymont et al. [1997] based on the time dependent MHD- 
code by Mikic et al. [1988]. The code uses a nonuniform 
mesh and the region of interested is embedded in a large 
computational domain to reduce the influence of the lat- 
eral boundaries. The method has been applied to extrap- 
olate the magnetic field above au active region by Jiao 
et al. [1997]. The computations have been carried out 
with a resolution of the order of 100'^ points. A super- 
computer was required for these computations that time 
(10 years ago), but due to the rapid increase of computer 
speed and memory within the last decade this restriction 
is very probably not valid anymore. 

Roumeliotis [1996] developed the so called stress and 
relax method. In this approach the initial potential field 
becomes disturbed by the observed transverse field com- 
ponent on the photosphere. The boundary conditions 
arc replaced in subsequently in several small steps and 
always relaxed with a similar MHD-relaxation scheme as 
described above towards a force-free equilibrium. The 
code by Roumeliotis [1996] has implemented a function 
w{x,y) which allows to give a lower weight to regions 
where the transverse photospheric field has been mea- 
sured with lower accuracy. Additional to the iterative 
equations as discussed above, the method includes a re- 
sistivity 77 (or diff'usivity) by adding a term r/j on the 
right hand site of Ohms law (22). This relaxes some- 
what the topological constrains of ideal MHD relaxation, 
because a finite resistivity allows a kind of artificial rccon- 
nection and corresponding changes of the initial potential 
field topology. The method has been tested with a force- 
free equilibrium found by Klimchuk and Sturrock [1992] 
and applied to an active region measured with the MSEC 
vector- magnctograph. 

The stress and relax method has been revisited by Val- 
ori et al. [2005]. Different from the earlier implementa- 
tion by Roumeliotis [1996] the new implementation uses 
directly the magnetic field, rather than the vector poten- 
tial in order to keep errors from taking numerical devia- 
tions from noisy magnetograms minimal. The solenoidal 
condition is controlled by a diffusive approach by Dedner 
et al. [2002] which removes effectively a numerically cre- 
ated finite divergence of the relaxed magnetic field. The 
new implementation uses a single stress step, rather than 
the multiple small stress used by Roumeliotis [1996] to 
speed up the computation. 

The single step stress and relax method is connected 
with a suitable control of artificial plasma flows by the 
Courant criterium. The authors reported that a multi- 
step and single-step implementation do not reveal sig- 
nificant differences. The numerical implementation is 
based on the time-depended fuU MHD-code 'AMRVAC 
by Keppens et al. [2003]. Valori et al. [2005] tested their 
nonlinear force-free implementation with a numerically 
constructed nonlinear force-free twisted loop computed 
by Tordk and Kliem [2003]. 
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2.4. Optimization approach 

The optimization approach has been developed in 
Wheatland et al. [2000]. The solution is found by mini- 
mizing the functional 



L= i [B-^|(V X B) X B|^ + |V-B|'] <fv. (29) 
Jv 

Obviously, L is bound from below by 0. This bound 
is attained if the magnetic field satisfies the force-free 

equations (l)-(3). 

By taking the functional derivatives with respect to 
some iteration parameter t we get: 

with 



F = V X 



+ 



[(V X B) X B] X B 



-V X 



52 

((V-B) B) X B 

52 



-n X (V X B) - V(f2-B) 
fi(V • B) + B} 



(31) 



n = B"^ [(V X B) X B - (V ■ B) B] (32) 

The surface term vanishes if the magnetic field vector 
is kept constant on the surface, e.g. prescribed from pho- 
tosphcric measurements. In this case L decreases mono- 
tonically if the magnetic field is iterated by 



f^.F^ (33) 

Let us remark that Fmhs as defined in Equation (27) 
and used for MHD-rclaxation is identical with the first 
term on the right-hand-side of Eq. (31), but Eq. (31) 
contains additional terms. 

For this method the vector field B is not necessarily so- 
lenoidal during the computation, but will be divergence- 
free if the optimal state with L = is reached. A dis- 
advantage of the method is that it cannot be guaranteed 
that this optimal state is indeed reached for a given ini- 
tial field and boundary conditions. If this is not the case 
then the resulting B will either be not force-free or not 
solenoidal or both. 

McTiernan has implemented the optimization ap- 
proach basically as described in Wheatland et al. [2000] 
in IDL (see Schrijver et al. [2006] for a brief descrip- 
tion of the McTiernan implementation.) This code al- 
lows the use of a non-uniform computational grid. In a 
code inter-comparison by Schrijver et al. [2006] the IDL 
optimization code by McTiernan was about a factor of 50 
slower compared to an implementation in parallelized C 
by Wiegelmann [2004] . To our knowledge McTiernan has 
translated his IDL-code into FORTRAN in the meantime 
for faster computation (personal communication on the 
NLFFF-workshop Palo- Alto, June 2006. See also Metcalf 
et al. [2007].) 

Several tests have been performed with the optimizar- 



WIEGELMANN: CORONAL MAGNETIC FIELDS. 



X- 13 



tion approach in Wiegelmann and Neuktrch [2003]. It 
has been investigated how the unknown lateral and top 
boundary influence the solution. The original optimiza- 
tion approach by Wheatland et al. [2000] has been ex- 
tended towards more flexible boundary-conditions, which 
allow ^ 7^ on the lateral and top boundaries. This has 
been made with the help of the surface integral term in 
(30) and led to an additional term ^ = fxG on the 
boundaries. This approached improved the performance 
of the code for cases, where only the bottom boundary 
was prescribed. No improvement was found for a slow 
multi-step replacement of the boundary and this possibil- 
ity has been abandoned in favour of a single stop method. 

It has been also investigated how noise influences the 
optimization code and this study revealed that noise in 
the vector magnetograms leads to less accurate nonlinear 
force- free fields. 

Wiegelmann [2004] has reformulated the optimization 
principle by introducing weighting functions One defines 
the functional 



L = I [wB-^|(V X B) X B|^ -h w |V -Bl^] d^a(34) 
Jv 

where w{x, y, z) is a weighting function. It is obvious 
that (for w > 0) the force- free equations (1-3) arc fulfilled 
when L is equal zero. Minimization of the functional (34) 
lead to: 



^=mF, (35) 

F = w F + (fia X B) X Vw (fib • B) Vw (36) 
f2a = B'^ [(V X B) X B] (37) 
fib = B"^ [(V-B)B], (38) 

with F as defined in (31). With w{x,y^z) = 1 this 
approach reduces to the Wheatland et al. [2000] method 
as described above. The weighting function is useful if 
only the bottom boundary data arc known. In this case 
we a buffer boundary of several grid points towards the 
lateral and top boundary of the computational box is in- 
troduced. The weighting function is chosen constant in 
the inner, physical domain and drop to with a cosine 
profile in the buffer boundary towards the lateral and top 
boundary of the computational box. In Schrijver et al. 
[2006] some tests have been made with different weight- 
ing functions for the force-free and solenoidal part of the 
functional (34), but the best results have been obtained 
if both terms got the same weight. The computational 
implementation involves the following steps. 

1. Compute start equilibrium (e.g. a potential field) 
in the computational box. 

2. Replace the bottom boundary with the vector mag- 
netogramm. 

3. Minimize the functional (34) with the help of Eq. 
(35). The continuous form of (35) guaranties a monoton- 
ically decreasing L. This is as well ensured in the dis- 
cretized form if the iteration step dt is sufficiently small. 
The code checks if L{t+dt) < L{t) after each time step. If 
the condition is not fulfilled, the iteration step is repeated 
with dt reduced by a factor of 2. After each successful 
iteration step we increase dt slowly by a factor of 1.01 to 
allow the time step to become as large as possible with 
respect to the stability condition. 
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4. The iteration stops if L becomes stationary. Sta- 
tionarity is assumed if ^/L < 1.0 • 10"'* for 100 consec- 
utive iteration steps. 

The program has been tested with the semi-analytic 
nonUnear force-free configuration by Low and Lou [1990] 
and Titov and Demoulin [1999] in Wiegelmann et al. 
[2006a]. The code has been appUed to extrapolate the 
coronal magnetic field in active regions in Wiegelmann 
et al. [2005b, a]. 

A finite clement optimization approach has been im- 
plemented by Inhester and Wiegelmann [2006] using the 
Whitney elements as for the Grad-Rubin code (which 
has been described above). The optimization method 
uses exactly the same staggered finite element grid as de- 
scribed above, which is different from the finite difference 
grids used in the earlier implementations by Wheatland 
et al. [2000]; Wiegelmann and Neukirch [2003]; Wiegel- 
mann [2004]. Another difference is that earlier imple- 
mentations discrctizcd the analytical derivative of the 
functional L (31), while the new code takes the numeri- 
cal more consistent derivative of the discrctizcd function 
L. All other implementations used a simple Landwcber 
scheme for updating the magnetic field, which is replaced 
here by an unpreconditioned conjugate gradient iteration, 
which at every time step performs an exact line search 
to the minimum of L in the current search direction and 
additional selects an improve search direction instead of 
the gradient of the functional L. To do so the Hessian 
matrix of the functional L is computed during every itera- 
tion step. An effective computation of the Hessian matrix 
is possible, because the reformulated function L{s) is a 
fourth order polynomial in B and all five polynomial co- 
efficients can be computed in one go. The code has been 
tested with Low and Lou [1990] and the result of twisted 
loop computations of the Grad-Rubin implementation on 
the same grid. 

The optimization code in the implementation of 
Wiegelmann [2004] has recently be extended towards us- 
ing a multi-scale implementation. The main difference 
from the original code are (see also Metcalf et al. [2007]): 

• The method is not full multigrid, but computes the 
solution on different grids only once, e.g., something like 
50^ 100^ 200^ 

• The main idea is to get a better (than potential field) 
start equilibrium on the full resolution box. 

• Solution of smaller grids are interpolated onto larger 
grids as initial state for the magnetic field in the compu- 
tational domain of the next larger box. 

The multiscale implementation has been tested as part 
of a codc-intcr-comparison test in Metcalf et al. [2007] 
with the help of solar-like reference model computed by 
van Ballegooijen [2004]; van Ballegooijen et al. [2007]. 

The optimization approach has recently be imple- 
mented in spherical geometry by Wiegelmann [2007] and 
tested with Low and Lou [1990]. The original longitudi- 
nal symmetric Low and Lou solution has been shifted by 
1/4 of a solar radius to test the code without any sym- 
metry with respect to the Suns surface. The numerical 
implementation is very similar as the Cartesian imple- 
mentation described in Wiegelmann [2004]. The spher- 
ical implementation converged fast for low latitude re- 
gions, but the computing time increased significantly if 
polar regions have been included. It has been suggested 
to implement the code on a so called 'Yin and Yang' grid 
as developed by Kageyama and Sato [2004] to reduce the 
computing time. The 'Yin and Yang' grid is suitable for 
massive parallelization, which is necessary for full sphere 
high resolution NLFFF-computations. 
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2.5. Boundary element or Greens function like 
method 

The boundary integral method has been developed by 
Yan and Sakurat [2000]. The method relates the mea- 
sured boundary values with the nonlinear force-free field 
in the entire volume by: 

where Ci — 1 for points in the volume and Ci = 1/2 for 
boundary points and Bo is the measured vector magnetic 
field on the photosphere. The auxiliary vector function 
is defined as 

V - Hiao- f ''O'^i^'^'^) cosjXyr) cos{X,r) \ 
- ^'^S 1^ Attv ' 47rr ' Attv J ^ 

and the \i, {i = x,y,z) are computed in the original 
approach by Yan and Sakurai [2000] with integrals over 
the whole volume, which define the Aj implicitly: 



Yi[XiBi-a^Bi-{VaxBi)]dV = (41) 

This volume integration, which has to be carried out 

for every point in the volume is certainly very time con- 
suming (a sixth order process). The \i have the same 
dimension as the magnetic field. The existence of the Ai 
has been confirmed for the senu-analytic field of Low and 
Lou [1990] by Li et al. [2004]. While the work of Li et al 
[2004] showed that one can find the auxiliary function Y 
for a given force-free field in 3D, the difficulty is that Y is 
a/-priori unknown if only the photospheric magnetic field 
vector is given. Yan and Sakurai [2000] proposed an iter- 
ative scheme to compute the auxiliary functions and the 
nonlinear force-free magnetic field selfconsistently. They 
use the approximate solution k on the right hand side of 
Eq. (39) to compute a better solution k + 1 by: 

where the initial guess for the magnetic field in the 
volume is B = and also the initial ^ = 0. In principle 
it would be also possible to compute a potential field first 
and derive the auxiliary functions for this field as done in 
Li et al. [2004] and iterate subsequently for the nonlin- 
ear force-free fields and the associated auxiliary functions 
with Eq. (42). This possibility has not been tried out to 
our knowledge until now, however. The method iterates 
the magnetic field until B and ^ converge. In an inter 
code comparison by Schrijver et al. [2006] one iteration 
step of (42) took about 80h for this method and only 
this one step was carried out without further iteration. 
This seems, however, not to be sufficient to derive an 
accurate nonlinear force-free solution. The method has 
been applied for the comparison with soft X-ray loops ob- 
served with YOHKOH by Wang et al. [2000]; Liu et al. 
[2002] and to model a magnetic flux robe by Yan et al. 
[2001a, b]. 

In a new implementation of the boundary element 
method by Yan and Li [2006] the auxiliary functions are 
computed itcratively with the help of a simplex method. 
This avoids the numerical expensive computation of the 
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volume integral (41). The boundary element method is 
still rather slow if a magnetic field has to be computed 
in an entire 3D domain. Different from other method, it 
allows, however, to evaluate the NLFFF-field at every ar- 
bitrary point within the domain from the boundary data, 
without the requirement to compute the field in an entire 
domain. This is in particular useful if one is interested 
to compute the NLFFF-field only along a given loop. 

He and Wang [2006] investigated the validity of the 
boundary integral representation for a spherical imple- 
mentation. The method has been tested with the lon- 
gitudinal invariant Low and Lou [1990] solution. The 
spherical implementation method of this method revealed 
reasonable results for smooth modestly nonlinear fields, 
but a poor convergence for complex magnetic field struc- 
tures and large values of a. 

3. How to deal with non-force-free 
boundaries and noise ? 

Given arbitrary boundary conditions of the magnetic 
field vector on the photosphere, the solution to the force- 
free equations in 3D may not exist. Nonlinear force-free 
coronal magnetic field models assume, however, that the 
solution exists. It is certainly possible and necessary to 
check after or during the computation if a solution has 
been found. In the following we will discuss what we can 
do if the measured photospheric data are incompatible 
with the assumption of a force-free coronal magnetic field. 

3.1. Consistency check of vector magnetograms 

We reexamine some necessary conditions with the pho- 
tospheric field (or bottom boundary of a computational 
box). These conditions have to be fulfilled in order to 
be suitable boundary conditions for a nonlinear force- 
free coronal magnetic field extrapolation. An a-priori as- 
sumption about the photospheric data is that the mag- 
netic flux from the photosphere is sufficiently distant from 
the lateral boundaries of the observational domain and 
the net flux is in balance, i.e.. 



Molodensky [1969, 1974]; Aly [1989]; Sakurai [1989] 
used the virial theorem to define which conditions a vec- 
tor magnetogram has to fulfill to be consistent with the 
assumption of a force-free field in the corona above the 
boundary. These conditions are: 

1. The total force on the boundary vanishes 




(43) 




(44) 



(45) 



2. The total 



torque on the boundary vanishes 




(46) 



(47) 



(48) 
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In an earlier review Aly [1989] has mentioned already 
that the magnetic field is probably not force-free in the 
photosphere, where B is measured because the plasma 
/3 in the photosphere is of the order of one and pres- 
sure and gravity forces are not negligible. The integral 
relations (44)- (48) are not satisfied in this case in the 
photosphere and the measured photospheric field is not 
a suitable boundary condition for a force-free extrapola- 
tion. Metcalf et al, [1995] concluded that the solar mag- 
netic field is not forco-frcc in the photosphere, but be- 
comes forco-frcc only at about 400fcm above the photo- 
sphere. Gary [2001] pointed out that care has to be taken 
when extrapolating the coronal magnetic field as a force- 
free field from photospheric measurements, because the 
force-free low corona is sandwiched between two regions 
(photosphere and higher corona) with a plasma ,9 w 1, 
where the force-free assumption might break down. An 
additional problem is that measurements of the photo- 
spheric magnetic vector field contain inconsistencies and 
noise. In particular the transverse components (say Bx 
and By) of current vector magnetographs include uncer- 
tainties. 

The force-free field in a domain requires the Maxwell 
stress (44)-(48) to sum to zero over the boundary. If 
these conditions are not fulfilled a force-free field can- 
not be found in the volume. A faithful algorithm should 

therefore have the capability of rejecting a prescription 
of the vector field at the boundary that fails to produce 
zero net Maxwell stress. A simple way to incorporate 
these conditions would be to evaluate the integrals (44)- 
(48) within or prior to the NLFFF-computation and to 
refuse the vector field if the conditions are not fulfilled 
with sufficient accuracy. Current codes do run, however, 
although if feeded with inconsistent boundary data, but 
they certainly cannot find a force-free solution in this 
case (because it docs not exist). This property of current 
codes docs, however, not challenge the trustworthiness 
of the algorithms, because the force-free and solcnoidal 
conditions are checked in 3D, e.g., with the help of the 
functional L as defined in (29). A non zero value of L 
(within numerical accuracy) tells the user that a force- 
free state has not been reached. In principle it would be 
possible that the codes do refuse to output the magnetic 
field in this case. For current codes this is not automati- 
cally controlled but responsibility of the user. 

Unfortunately current measurements of the magnetic 
field vector arc only available routinely in the pho- 
tosphere, where we have a finite (3 plasma and non- 
magnetic forces might become important. The force-free 
compatibility conditions (44)- (48) are not fulfilled in the 
photosphere, but they should be fulfilled in the low /J 
chromospheric and coronal plasma above. The question 
is if we still can use the photospheric measurements to 
find suitable consistent boundary conditions for a non- 
linear force-free modelling. Such an approach has been 
called preprocessing of vector magnetograms. 

3.2. Preprocessing 

The preprocessing routine has been developed by 
Wiegelmann et al. [2006b]. The integral relations (44)- 
(48) have been used to define a 2D functional of quadratic 
forms: 



Lprep = fJ-iLi + fj.2L2 + HsLs + iuLi (49) 

where 
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Li = 



B^B^j + ByB.j + Bl -Bl- BlJO 
{^x{Bl-Bl-Bl)^ +{^y[Bl-Bl-Bl)^ 



L4 = 



+ l^yB^B^ - xByB^j 

^ ^ {Bx — Bxobs)^ + ^ ^ {By — Byobs)'^ 
. P 

+ '^{B,- B^obsf 

p 

J2 i^Bxf + {AByf + {AB^f 



(51) 



(52) 



(53) 



The surface integrals are here replaced by a summa- 
tion over all grid nodes p of the bottom surface grid 
and the differentiation in the smoothing term is achieved 
by the usual 5-point stencil for the 2D-Laplace opera- 
tor. Each constraint L„ is weighted by a yet undeter- 
mined factor fin- The first term (n=l) corresponds to 
the force-balance conditions (44)-(45), the next (n=2) to 
the torque- free condition (46)- (48). The following term 
(n=3) ensures that the optimized boundary condition 
agrees with the measured photospheric data and the last 
terms (n=4) controls the smoothing. The 2D-Laplace 
operator is designated by A. The aim of the preprocess- 
ing procedure is to minimize Lprep so that all terms Ln 
if possible are made small simultaneously. A strategy 
on how to find the optimal yet undefined parameters /i„ 
is described in Wiegelmann et al. [2006b]. As result of 
the preprocessing we get a data set which is consistent 
with the assumption of a force-free magnetic field in the 
corona but also as close as possible to the measured data 
within the noise level. 



4. Code testing and code comparisons 

Newly developed codes for the extrapolation of nonlin- 
ear force-free fields from boundary data have to be tested 
before they are applied to measurements. In principle any 
analytical or numerically created solution of the force-free 
equations (l)-(3) can be used as a reference case. One 
cuts a plane (artificial photosphere, bottom boundary) 
out of the 3D reference solution ^ and uses the above de- 
scribed extrapolation codes to reconstruct the magnetic 
field. The result of this extrapolation is then compared 
with the reference to rate the quality of the reconstruc- 
tion. Unfortunately, it is very hard to find a truly non- 
linear 3D solution of (l)-(3) analyticallj' and very few 
solutions are known. Low and Lou [1990] (LL) found a 
class of solutions which have become a standard refer- 
ence for testing NLFFF-extrapolation codes. LL found 
axisymmetric equilibria which are separable in spherical 
coordinates. They are self-similar in the radial coordi- 
nate, and the polar angle dependence is determined from 
a nonlinear eigenvalue equation. The symmetry is bro- 
ken by cutting out a rectangular chunk of the solution by 
using a Cartesian coordinate system which is shifted and 
rotated with respect to the original coordinate system in 
which the LL equilibria are calculated. The parameters of 
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the LL solutions and the parameters of the new Cartesian 
coordinate system allow for a large number of different 
situations which can be used for tests. The original ax- 
isymmetric spherical LL solution has also been used (with 
and without symmetry breaking by shifting the origin 
of the coordinate system) to test spherical NLFFF pro- 
grams. To our knowledge all recent implementations of 
the described NLFFF approaches have been tested with 
LL, cither immediately in the original code-describing pa- 
pers or in subsequent works, e.g., in a blind-algorithm 
test within the NLFFF-consortium, as described below. 

The MHD-rolaxation method and the optimization ap- 
proach have been compared by Wiegelmann and Neukirch 
[2003]. Both methods have been applied to the Low and 
Lou [1990] equilibrium with exactly the same finite dif- 
ference grid. The iterative equations for MHD-relaxation 
and optimization have both the form 4^ = /iF but the 
structure of F is more complicated for optimization than 
for MHD-rolaxation. The MHD-rclaxation term is in- 
deed identical with the first term of the optimization 
approach. While MHD-relaxation minimizes only the 
Lorentz-force, the optimization does additional minimize 
V-B, while a decreasing magnetic field divergence during 
MHD-relaxation (as seen in Wiegelmann and Neukirch 
[2003]) is the result of numerical diffusion. Despite the 
numerical overhead in computing F for the optimization 
code, optimization provided more accurate results and 
faster convergence. 

A practical advantage of the MHD-approach is that 
several time-dependent MHD-codes are well known and 
established and can be used for the force-free relaxation 
discussed here. The inclusion of non-magnetic forces 
like pressure gradients and gravity looks straight for- 
ward for the MHD approach. Other methods are usu- 
ally developed with the only task of computing nonlin- 
ear forcc-froo coronal magnetic fields, also a generaliza- 
tion towards magnctohydrostatic and stationary MHD- 
equilibria is possible and has been done for the opti- 
mization approach [see Wiegelmann and Inhester, 2003; 
Wiegelmann and Neukirch, 2006]. Another advantage of 
using time-dependent MHD-codes for rela^cation is that 
the computed force-free equilibrium can be used on the 
same grid and with the same code as initial state for 
time-dependent MHD-simulations. One can, in princi- 
ple, use the force-free equilibria computed with any of 
the described method as initial state for time-dependent 
MHD-simulation, but having the initial equilibrium state 
already directly on the MHD-grid might be very handy, 
because no further adjustments are needed. 

4.1. The NLFFF-consortium 

Since the year 2004 activities arc ongoing to bring 
NLFFF-modelers together and to compare the different 
existing codes. A workshop series has been organized 
for this aim by Karel Schrijver and three workshops took 
place so far from 2004-2006. The next workshop is planed 
for June 2007. As we have been asked to summarize 
the workshop results on the CSWM-meeting, we give 
also a very brief overview in the corresponding special 
issue-paper here. The main results of the first two work- 
shops have been published in Schrvjver et al. [2006]. In 
this paper six difl^erent NLFFF-implcmcntations (Grad- 
Rubin codes of Amari et al. [1999] and Wheatland [2004], 
MHD-relaxation code of Valori et al. [2005] , optimization 
codes by McTiernan and Wiegelmann [2004], boundary 
element method by Yan and Sakurai [2000]) have been 
compared. The codes have been tested in a blind algo- 
rithm test with the help of the semi-analytic equilibrium 
by Low and Lou [1990] in two cases. In case I all six 
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boundaries of a computational box have been described 
and in case II only the bottom boundary. The compar- 
ison of the extrapolation results with the reference solu- 
tion has been done qualitatively by magnetic field line 
plots (Shown here in Fig. 3 for the central region of case 
II) and quantitatively by a number of sophisticated com- 
parison metrices. All NLFFF-fields agreed best with the 
reference field for the low lying central magnetic field 
region, where the magnetic field and electric currents 
are strongest and the influence of the boundaries low- 
est. The code converged with speeds that differed by a 
factor of one million per iteration steps The fastest- 
converging and best-performing code was the Wheatland 
et al. [2000] optimization code as implemented by Wiegel- 
mann [2004] . Recent implementations of the Grad- Rubin 
code by Amari et al. [2006]; Inhester and Wiegelmann 
[2006] and a new implementation of the upward integra- 
tion method by Song et al. [2006] did not participate in 
the blind-algorithm intcr-comparison by Schrijver et al. 
[2006] , but these three new codes have been tested by the 
authors with similar measures and revealed similar accu- 
racy as the best performing codes in the blind algorithm 
test. It seems that the somewhat more flexible boundary 
conditions used in the Grad-Rubin approaches of Amari 
et al. [2006] and Inhester and Wiegelmann [2006] are re- 
sponsible for the better performance compared to the 
earlier implementation by Amari et al. [1999], which has 
been used in the blind algorithm test. 

The widely used LL-cquilibrium contains a very 
smooth photospheric magnetic field and an extended cur- 
rent distribution. It is therefore also desirable to test 
NLFFF-codes also with other, more challenging bound- 
ary fields, which are less smooth, have localized current 
distribution and to investigate also the effects of noise 
and effects from non force-free boundaries. A somewhat 
more challenging reference case is the equilibrium found 
by Titov and Demoulm [1999] (TD). Similar as LL, the 
TD equilibrium is an axisymmetric equilibrium. The TD- 
model contains a potential field which is disturbed by a 
toroidal nonlinear force-free current. This equilibrium 
has been used for testing the MHD-rela:xation code (Val- 
ori and Kliem, private communication) and the optimizar 
tion code in Wiegelmann et al. [2006a]. 

Any numerically created NLFFF-model might be suit- 
able for code testing, too. It is in particular inter- 
esting to use models, which arc partly related on ob- 
servational data. Very recently van Ballegooijen et al. 
[2007] used line-of-sight photospheric measurements from 
SOHO/MDI to compute a potential field, which was then 
disturbed by inserting a twisted flux robe and relajced 
towards a nonlinear force-free state with a magnetofric- 
tional method as described in van Ballegooijen [2004]. 
The van Ballegooijen et al. [2007] model is not force- 
free in the entire computational domain, but only above 
a certain height above the bottom boundary (artificial 
chromosphere). On the lowest boundary (photosphere) 
the model contains significant non-magnetic forces. Both 
the cliromospheric as well as the photospheric magnetic 
field vector from the van Ballegooijen et al. [2007]-modcl 
have been used to test four of the recently developed ex- 
trapolation codes (One Grad-Rubin method, one MHD- 
relaxation code and two optimization approaches) in a 
second blind algorithm test by Metcalf et al. [2007]. 
While the I. NLFFF-consortium paper {Schrijver et al. 
[2006]) used a domain of just 64'^ pixel, the 11. paper used 
a computational domain of 320 x 320 x 258 pixel and 
modern NLFFF-codes where able to compute the non- 
linear force-free field in such relatively large boxes within 
a few hours for a moderate parallelization on only 1-4 
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processors and a memory requirement of 2.5 — 4 GB of 
Ram. This very recent code-comparison shows a major 
improvement regarding computing time and suitable grid 
sizes within less than three years. On the first NLFFF- 
consortium meeting in 2004, box sizes of some 64^ have 
been a kind of standard or computing times of some two 
weeks have been reported for 150^ boxes. We briefly sum- 
marize the results of Metcalf et al. [2007] as: 

• NLFFF-oxtrapolations from chromosphcric data re- 
cover the original reference field with high accuracy. 

• When the extrapolations are applied to the photo- 
spheric data, the reference field is not well recovered. 

• Preprocessing of the photospheric data improve the 
result, but the accuracy is still lower as for extrapolations 
from the chromosphere. 



5. Conclusions and Outlook 

Within the last few years the scientific community 
showed a growing interest into coronal magnetic fields. 
^ The development of new ground based and space born 
vector magnetographs provide us measurements of the 
magnetic field vector on the suns photosphere. Accompa- 
nied from these hardware development, software has been 
developed to extrapolate the photospheric measurements 
into the corona. Special attention has recently been 
given to nonlinear force-free codes. Five different numer- 
ical approaches (Grad-Rubin, upward integration, MHD- 
rela;xation, optimization, boundary elements) have been 
developed for this aim. It is remarkable that new codes 
or major updates of existing codes have been published 
for all five methods within the last two years, mainly 
in the last year (2006). A workshop series (NLFFF- 
consortium) since 2004 on nonlinear force-free fields has 
recently released synergy effects, by bringing modelers of 
the different numerical implementations together to com- 
pare, evaluate and improve the programs. Several of the 
most recent new codes and utility programs (e.g. pre- 
processing) have at least been parth^ inspired by these 
workshops. The new implementations have been tested 
with the smooth semi-analytic Low-Lou-cquilibrium and 
showed reasonable agreement with this reference field. 
While all methods aim for a reconstruction of the coro- 
nal magnetic field from the photospheric magnetic field 
vector, the way how these measurements are used to pre- 
scribe the boundaries of the codes is different. 

• MHD-relaxation and optimization use Bxo, Byo, B^o 
on the bottom boundary. This over-determines the 
boundary value problem. Both methods axe closely re- 
lated and compute the magnetic field in a computational 
box with 



f=.F, ,54) 

where the structure of F is somewhat different (the op- 
timization approach has more terms) for both meth- 
ods. Usually a potential field is used as initial state for 
both approaches, also the use of a linear force-free ini- 
tial state is possible. Recently a multiscale version of 
optimization has been installed, which uses a low reso- 
lution NLFFF-field as input for higher resolution com- 
putations. Specifying the entire magnetic field vector on 
the bottom boundary is an over-imposed problem and a 
unique NLFFF-field (or a solution at all) requires that the 
boundary data fulfill certain consistency criteria. A re- 
cently developed preprocessing-routine helps to find suit- 
able consistent boundary data from inconsistent photo- 
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spheric measurements. Earlier and current comparisons 
showed a somewhat higher accuracy for the optimization 
approach. A practical advantage of the MHD-approach 
is that in principle any available time-dependent MHD- 
code can be adjusted to compute the NLFFF-field. 

• The Grad-Rubin approach uses B^o and the distri- 
bution of a computed with Eq. (8) for one polarity, 
which corresponds to well posed mathematical problem. 
A practical problem is that the computation of a requires 
numerical differences of the noisy and forced transverse 
photospheric field Bxo, Byo with (7) leading to inaccura- 
cies in the normal electric current distribution and in a. 
For smooth semi-analytic test cases this is certainly not a 
problem, but real data require special attention (smooth- 
ing, preprocessing, limiting a 7^ to regions where B^o is 
above a certain limit) to derive a meaningful distribution 
of a. While the method requires only a for one polarity, 
the computation from photospheric data provide a for 
both polarities. We are not aware of any tests on how 
well NLFFF-solutions computed from a prescribed on 
the positive and negative polarity coincide. It is also un- 
clear how well the computed transverse field components 
on the bottom boundary agree with the measured values 
^ of Bxo , Byo . More tests on this topics are necessary, 
including the recently installed possibility to prescribe a 
for both polarities and adjust the boundary by a weighed 
average of a on both polarities to fulfill Eq. (6). As ini- 
tial state the Grad-Rubin method uses a potential field, 
which is also true for MHD-relaJcation and optimization. 

• The upward integration and the boundary element 
method prescribe both all components of the bottom 

boundary magnetic field vector and the a distribution 
computed with Eq. (8). This approach over imposes 
the boundary and Bxo, Byo, Bzo and a have to be con- 
sistent which each other and the force-free assumption. 
This is certainly not a problem at all for smooth semi- 
analytic test equilibria and strategies to derive consistent 
boundary data from measured data have been developed 
recently. Different from the three approaches discussed 
above, upward integration and boundary element meth- 
ods do not require to compute first an initial potential 
field in the computational domain. It is well known that 
the upward integration method is based on an ill-posed 
problem and the method has not been considered for sev- 
eral years, but a recent implementation with smooth ana- 
lytic functions might help to regularize this method. First 
tests showed a reasonable results for computations with 
the smooth semi-analytic Low-Lou solution. 
The boundary element method has the problem to be 
very slow and an earlier implementation of this method 
could not reach a converged state for a 64'^ boxed used in 
the I. NLFFF-consortium paper due to this problem. A 
now 'direct boundary method' has been developed, which 
seems to be faster than the original 'boundary clement 
method', but still slower compared with the four other 
NLFFF-approaches if the task is to compute a 3D mag- 
netic field in an entire 3D-domain. Different from all 
other described methods the boundary element approach 
allows to compute the nonlinear force-free field vector 
at any arbitrary point above the boundary and it is not 
necessary to compute the entire 3D-field above the pho- 
tosphere. This might be a very useful feature if one is 
interested in computing the magnetic field only along a 
single loop and not interested in an entire active region. 
The new implementations of upward integration and 
boundary element method show both reasonable results 
for first tests with the smooth semi-analytic Low and Lou 
equilibrium. Further tests with more sophisticated equi- 
libria, e.g. a solar-like test case as used in the II. NLFFF- 
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consortium paper would be useful to come to more sound 
conclusions regarding the feasibility of these methods. 

Most of the efforts done in nonlinear force-free mod- 
elling until now concentrated mainly on developing these 
models and testing their accuracy and speed with the 
help of well known test configuration. Not too many ap- 
plications of nonlinear force-free models to real data are 
currently available, from which we learned new physics. 
One reason was the insufficient access to high accuracy 
photosphcric vector magrictograms and a second one were 
limitations of the models. Force-free field extrapolation 
is a mere tool, if properly employed on vector magne- 
tograms, it can help to understand physical, magnetic 
field dominated processes in the corona. Both the com- 
putational methods as well as the accuracy of required 
measurements (e.g. with Hinode, SDO) are rapidly im- 
proving. Within the NLFFF-consortium we just started 
(since april 2007) to apply the different codes to compute 
nonlinear force-free coronal magnetic fields from Hinode 
vectormagnetograms. This project might provide us al- 
ready some new insights about coronal physics. 

To conclude, wc can say that the capability of Carte- 
sian nonlinear force-free extrapolation codes has rapidly 
increased in recent years. Only three years ago most 
codes run usually on grids of about 64^ pixel. Recently 
developed or updated codes (Grad-Rubin by Wheatland, 
MHD-rclaxation by Valori, optimization by Wicgclmann, 
optimization by McTicrnan) have been applied to grids 
of about 300^ pixel. Although this increase of trace- 
able grid sizes is certainly encouraging, the resolution of 
current and near future vector magnetographs (which of 
course measure only data in 2D!) is significantly higher. 
We should keep in mind, however, that the currently 
implemented NLFFF-codes have been only moderately 
parallelized using only a few processors. The CSWM- 
conference, where this paper has been presented, took 
place at the 'Earth simulator' in Yokohama, which con- 
tains several thousands of processors used for Earth- 
science computer simulations. An installation of NLFFF- 
codes on such massive parallel computers (which has been 
briefiy addressed on NLFFF-consortium meetings) com- 
bined with adaptive mesh refinements might enable dras- 
tically improved grid sizes. One shoxild not underesti- 
mate the time and effort necessary to program and in- 
stall such massive parallelized versions of existing codes. 
As full disk vectormagnetograms will become available 
soon (SOLIS, SDO/HMI) it is also an important task to 
take a spherical geometry into account. First steps in 
this direction have been carried out with the optimiza- 
tion and boundary element methods. Spherical NLFFF- 
geometries are currently still in it's infancy and have been 
tested until now only with smooth semi-analytic Low and 
Lou equilibria and require further developments. 

Attention has also recently been drawn to the prob- 
lem that the coronal magnetic field is force-free, but the 
photosphcric one is not. Tests with extrapolations from 
solar-like artificial photosphcric and chromospheric mea- 
surements within the II. NLFFF-consortium paper re- 
vealed that extrapolations from the (force-free) chromo- 
spheric field provide significantly better results as extrap- 
olations using directly the (forced) photosphcric field. 
Applying a pre-processing program on the photospheric 
data, which effectively removes the non-magnetic forces, 
leads to significantly better results, but they are not as 
good as by using the chromospheric magnetic field vec- 
tor as boundary condition. An area of current research 
is the possibility to use chromospheric images to improve 
the preprocessing of photospheric magnetic field mesisure- 
ments. Improvements in measuring the chromospheric 
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magnetic field directly [e.g. Lagg et al, 2004] might fur- 
ther improve to find suitable boundary conditions for 
NLFFF-extrapolations. Force- free extrapolations are not 
suitable, however, to understand the details of physical 
processes on how the magnetic field evolves from the 
forced photosphere into the chromosphere, because non- 
magnetic forces are important in the photosphere. For a 
better understanding of these phenomena more sophisti- 
cated models which take pressure gradients and gravity 
(and maybe also plasma flow) into account arc required. 
Some first steps have been done with a generalization of 
the optimization method by Wiegelmann and Neukirch 
[2006], but such approaches are still in their infancy and 
have been tested so far only with smooth MHD-equilibria. 
It is also not entirely clear how well necessary information 
regarding the plasma (density, pressure, temperature, 
flow) can be derived from measurements. Non-magnetic 
forces become important also in quiet sun regions (Schri- 
jver and van Ballegooijen [2005]) and in the higher lay- 
ers of the corona, where the plasma /? is of the order of 
unity. Coronagraph measurements, preferably from two 
viewpoints as provided by the STEREO-mission, com- 
bined with a tomographic inversion might help here to 
get insights in the required 3D structure of the plasma 
density. One should also pay attention to the combi- 
nation of extrapolation methods, as described here, with 
measurements of the Hanle and Zecman effects in coronal 
lines which allows the reconstruction of the coronal mag- 
netic field as proposed in feasibility studies of vector to- 
mography by Kramar et al. [2006] ; Kramar and Inhester 
[2006]. Other measurements of coronal features, e.g., 
coronal plasma images from two STEREO-viewpoints, 
can be used for observational tests of coronal magnetic 
field models. Using two viewpoints provide a much more 
restrictive test of models as images from only one view 
direction. While a nonlinear force-free coronal magnetic 
field model helps us to derive the topology, magnetic field 
and electric current strength in coronal loops, they do not 
provide plasma parameters. One way to get insights re- 
garding the coronal plasma is the use of scaling laws to 
model the plasma along the reconstructed 3D field lines 
and compare correspondent artificial plasma images with 
real coronal images. Schrijver et al. [2004] applied such 
an approach to global potential coronal magnetic fields 
and compared simulated and real coronal images from 
one viewpoint. A generalization of such methods towards 
the use of more sophisticated magnetic field models and 
coronal images from two STEREO-viewpoints will prob- 
ably provide many insights regarding the structure and 
physics of the coronal plasma. An important challenge is 
for example the coronal heating problem. The dominat- 
ing coronal magnetic field is assumed to play an impor- 
tant role here, because magnetic field configuration con- 
taining free energy can under certain circumstances re- 
connect {Priest [1996, 1999]) and supply energy for coro- 
nal heating. Priest et al. [2005] pointed out that mag- 
netic reconnection at separators and separatrices plays 
an important role for coronal heating. Nonlinear force- 
free models can help here to identify the magnetic field 
topology, magnetic null points, separatrices and localized 
strong current concentration. While magnetic reconnec- 
tion [see e.g. Priest and Schrijver, 1999] is a dynamical 
phenomenon, the static magnetic field models discussed 
here can help to identify the locations favourable for re- 
connection. Time sequences of nonlinear force-free mod- 
els computed from corresponding vector magnetograms 
will also tell wether the topology of the coronal magnetic 
field has changed due to reconnection, even if the physics 
of reconnection is not described by force-free models. So- 
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phisticated 3D coronal magnetic field models and plasma 
images from two viewpoints might help to constrain the 
coronal heating function further, which has been done so 
far with plasma images from one viewpoint [Aschwanden, 
2001a, b, by using data from Yokoh, Soho and Trace]. 
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Notes 

1. The authors used a somewhat different nomenclature: The 
upward integration method was called 'progressive ex- 
tension method' and the Grad-Rubin method 'iterative 
method'. That time the term 'iterative method' was reason- 
able because Grad-Rubin was the only iterative approach 
available, but now, 17 years later, several other iterative 
methods are available to compute nonlinear force-free fields. 

2. As we will sec below the 'optimization' approach leads to 
a similar iteration equations for the magnetic field, but a 
different artificial driving force F. 

3. For the pure task of code testing it is also acceptable to use 
all six boundaries of the reference solution. These kind of 
data are not available for real solar cases of course. 

4. The codes run on different machines, have been written in 
different programming languages and used different compil- 
ers. A real test of the exact computing time would comprise 
a proper operation count, e.g., the number of fixed point 
additions and multiplications per iteration step. 

5. Publications containing the phrase 'coronal magnetic fields' 
in title or abstract have been cited less than about 50 times 
per year until the early 1990th and this number increased 
to about 150 citations per year in 2004. A peak year was 
2006 (last year) with more than 300 citations. (Source: ISI 
Web of Knowledge, march 2007) 

6. In principle BxQ,Bya may have an additional field 
(BxO, Byo) + {dxdy)4> without making a difference for a 
and hence for the Grad-Rubin result. 
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Figure 1. Linear force- free field model for NOAA 7986 
with the best fit for a. The upper panel shows a group of 
loops with Q = 2.5 and the lower panel another group of 
loops with a = —2.0. The different optimal values of the 
linear force-free parameter within one active region are a 
contradiction to the linear assumption (a constant) and 
tell us that a consistent modelling of this active regions 
requires a nonlinear force-free approach. [This figure has 
been originally published as Fig. 7 in Wiegelmann and 
Neukirch [2002]. Used with permission of Springer.] 



X - 32 



WIEGELMANN: CORONAL MAGNETIC FIELDS. 




Measured loops Potential field 




Linear Force -Free Nonlinear Force- Fre 



Figure 2. Magnetic field structure of tlie newly devel- 
oped active region NOAA 9451. Direct measurements 
of the field have been compared with a potential, linear 
and nonlinear force-free model. Best agreement has been 
found for the nonlinear model. [This figure has been orig- 
inally published as part of Fig. 1 in Wiegelmann et al. 
[2005b]. Used with permission of A & A.] 
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Figure 3. Evaluation of six nonlinear forcc-frcc codes. 
The reference solution (a) has been compared with ex- 
trapolations with 'optimization' (b,c); 'MHD-relaxation' 
(d), 'Grad-Rubin' (e,f), 'Boundary element' (g). For 
comparison a linear-force-free (h) and potential field (i) 
are shown, too. The images show the central domain 
of the model. Only the bottom boundary has been used 
provided for the extrapolation. [This figure has been orig- 
inally published as Fig. 4 in Schrijver et al. [2006]. Used 
with permission of Springer.] 



